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We perform numerical studies of the BFKL and CCFM equations for the unintegrated gluon distribution 
supplemented with an absorptive boundary which mimics saturation. For the BFKL equation, this pro- 
cedure yields the same results for the saturation momentum and the gluon distribution above saturation 
as the non-linear BK equation, for both fixed and running coupling, and for all the considered energies. 
This similarity goes beyond expectations based on the correspondence with statistical physics, which 
i-P ' hold only for fixed coupling and asymptotically high energies. For the CCFM equation, whose non-linear 

generalization is not known, our method provides the first study of the approach towards saturation. We 
find that, in the running-coupling case, the CCFM and BFKL predictions for the energy dependence of 
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the saturation momentum are identical within our numerical accuracy. A similar saturation boundary 
could be easily implemented in the CCFM-based Monte Carlo event generators, so like CASCADE. 



1. Introduction 
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qq , The imminent high-energy experiments at LHC will considerably enlarge the phase-space where the 

unitarity corrections to QCD interactions, like gluon saturation and multiple scattering, are expected to 
be important. Such corrections should in particular influence some 'hard' observables, like jet production 
at forward rapidities, whose theoretical description lies within the realm of perturbative QCD. The jets 
0^ ! to be measured at LHC will carry relatively large transverse momenta Q > 10 GeV, but because of 

the high-energy kinematics, their description may go beyond the standard pQCD formalism at high Q 2 
- the DGLAP evolution |l| of the parton distributions together with the collinear factorization of the 
hadronic cross-sections. Rather, the high-energy evolution, of the BFKL or CCFM Q type, and the 
associated fc^-factorization scheme should prevail whenever the energy logarithms In s ~ Y are larger 
than the momentum ones, InQ 2 . Besides, this evolution is expected to be amended by non-linear effects 
reflecting gluon saturation and multiple scattering, whose theoretical description within pQCD has been 
given only recently [j, 0, 0, [H and (in full rigor) only in the leading-order approximation (see Refs. 
[3, E3] for recent extensions to running coupling) . Such effects can make themselves felt even at relatively 
large momenta Q, well above the saturation momentum Q s (the characteristic scale for the onset of 
unitarity corrections), via phenomena like geometric scaling 11, E2> EU EH EH, ES E3], which reflect 



the change in the unintegrated gluon distribution at high k± » Q s due to saturation at low k± < Q s . 
The saturation scale Q s grows, roughly, like a power of the energy: Q 2 ~ s A with A ~ 0.25 from fits to 
HERA data [11,0, 

which are also supported by next -to-leading order theoretical calculations [14j . For 



forward jet production in proton-proton collisions at LHC, Q s is expected in the ballpark of 2 to 3 GeV. 
Besides, much higher values of Q s can be effectively reached by focusing on 'hot spots' (partons at 
large x and high Q 2 , which develop their own gluon cascades and saturation momentum) in some rare 
events, so like Mueller-Navelet jets. 

In view of the above, it becomes important and urgent to provide realistic, quantitative, predictions 
for the effects of saturation on relatively hard (Q 2 ^> Q 2 ) observables at LHC. Besides making the whole 
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perturbative approach fully justified, the restriction to relatively hard momenta entails some important 
simplifications, which are essential for the strategy that we shall propose in this Letter. 

First, this implies that one can neglect the complex many-body correlations which develop at satu- 
ration (Q < Qs), but focus on the gluon phase-space density (or 'unintegratcd gluon distribution') alone. 
Hence, the standard /cy-factorization of the hadronic cross-sections still applies, but with modified un- 
integrated gluon distributions, which reflect saturation. This opens the possibility to include the effects 



of saturation within Monte-Carlo event generators based on fc^-factorization, so like CASCADE [21|, 
which relies on the CCFM evolution. 

Second, this means that the saturation effects in the gluon distribution and, in particular, the 
saturation momentum itself can be computed without a detailed knowledge of the non-linear dynamics 
responsible for unitarization. Rather, they are fully determined by the linear evolution if the latter 
is supplemented with an absorptive boundary condition at low momenta, whose position is energy- 
dependent (as this mimics the saturation momentum) and is self-consistently determined by the evolution 



12 . Il3l . Il4j . This property is both important and highly non-trivial. It is important as it allows one to 



perform studies of saturation even for evolution equations whose non-linear generalizations are not known, 
so like the CCFM evolution, and also the BFKL evolution beyond the leading-order approximation. It 
is moreover non-trivial since the high- energy evolution is non-local in transverse momenta, hence the 
growth in the gluon distribution at high momenta k± 3> Q s could be well feed by radiation from lower 
momenta k± < Q s . This is clearly the case for the linear evolution with running coupling, in which the 
gluon distribution grows faster in the infrared (where the coupling is stronger) and then acts as a source 
for radiating gluons at higher momenta. 

Yet, at least for a fixed coupling and for asymptotically high energies, it has been firmly established 
that the high-energy evolution towards saturation is driven by the linear evolution in the dilute tail of 
the gluon distribution at high transverse momenta k± 3> Q s - The respective argument is based on a 
correspondence between high-energy QCD and statistical physics 15|, [22J, which however fails to apply 
in the more realistic case of a running coupling [231 ]. For that case, our subsequent numerical results will 
provide the first unambiguous evidence that the evolution of the saturation front in QCD is indeed driven 
by the linear part of the evolution, at all energies, including the preasymptotic ones. Most precisely, we 
shall find that the BFKL equation with saturation boundary provides exactly the same results for the 
saturation momentum Q S (Y) and for the gluon distribution at k± > Q S (Y) as the non-linear Balitsky- 
Kovchegov (BK) equation 0, Q , for both fixed and running coupling, and for all the considered rapidities 
Y < 120 — including lower values Y < 14, as relevant for the phenomenology at LHC. 

But the BFKL equation and its non-linear, BK, extension, to be discussed in Sect. 2, will merely 
serve as a playground to test our method for numerically implementing the saturation boundary condition 
within a generic linear evolution. Our main interest is rather in the CCFM evolution, that we shall discuss 
in Sect. 3, and which stays at the basis of Monte-Carlo event generators [2l|. There are at least two reasons 
why the CCFM evolution is a privileged tool in that respect. First, it takes into account the quantum 
coherence between successive emissions, leading to angular ordering in the parton cascades. This allows 
for a more realistic description of the particle distribution in the final state as compared to the BFKL 
evolution. (The latter only guarantees the correct treatment of inclusive quantities like the unintegrated 
gluon distribution.) Second, the CCFM equation provides an interpolation between BFKL dynamics at 
small x and (approximate) DGLAP dynamics at larger x. Of course, when discussing saturation we are a 
priori interested in the small-x region, and there is little doubt that, for asymptotically high energies, the 
approach towards saturation is more correctly described by the BFKL formalism. However, the energies 
to be available at LHC are far from asymptotia, and a better treatment of the transverse momentum 
ordering within the parton cascades, as implicit in the CCFM formalism, is probably essential for most 
of the "small-x" phase-space to be experimentally accessible at LHC. 

In this Letter, we shall limit ourselves to applications of the CCFM and BFKL evolutions of the 
unintegrated gluon distribution, so it will be meaningful to compare their results. To achieve a faster 
numerical convergence, we shall solve the respective differential equations, rather than building Monte- 
Carlo codes. But our prescription for implementing the saturation boundary condition is straightforward 
to apply within a Monte-Carlo event generator, so like CASCADE. Our present analysis can be viewed 
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as a first step towards building an event generator which includes saturation |24[. The CCFM equation 
is considerably more involved than the BFKL one (in particular, its solution involves one additional 
variable: a maximum angle) , so in order to be able to follow this evolution up to relatively high rapidities 
we shall consider a slightly simplified version of it. The simplification is obtained by rewriting the CCFM 
equation in a 'more inclusive form', i.e., by using the virtual, Sudakov-like, terms to cancel some of 
the real gluon emissions [25|, [26|. This procedure also involves some kinematical approximations, which 
however are in the spirit of the CCFM formalisuQ. By solving the ensuing equation with saturation 
boundary condition, we shall for the first time study the onset of unitarity corrections within the CCFM 
evolution. One of our most interesting results is that, in the running coupling case, the CCFM evolution 
in the presence of saturation provides almost identical results (for the saturation momentum and the 
gluon spectrum above Q s ) as the respective BFKL evolution. 

2. BFKL evolution with absorptive boundary 

In this section we shall explain our method for effectively implementing saturation within a unitarity- 
violating linear evolution on the example of the BFKL equation 0]. This is interesting since the corre- 
sponding non-linear equation which obeys unitarity is known as well, the BK equation 0, [B[ , and thus 
it can be used to test our method. Although all our numerical studies will be performed in (transverse) 
momentum space (referring to the unintegrated gluon distribution), it is more convenient to explain our 
method in coordinate space. Then, the BK equation describes the high-energy evolution of the scattering 
amplitude T(Y, r) of a small quark-antiquark dipole with transverse size r. We shall assume the target to 
be infinite and homogeneous in transverse directions, so we can ignore the impact-parameter dependence 
of the scattering amplitude and average over angles. The corresponding equation reads 



dY 2tt J z 2 (r-z) 2 

{ - T(Y, r) + T(Y, z) + T(Y, \r - z\) - T{Y, z)T{Y, \r-z\)). 

Here a s = a s N c /n : and z and r — z are the transverse sizes of the two dipoles into which the parent 
dipole r has dissociated before scattering off the target. The last term, quadratic in T, in the r.h.s. of 
the equation describes multiple scattering (the simultaneous scattering of both daughter dipoles) and is 
responsible for unitarization. 

With this last term omitted, (|2.ip reduces to the BFKL equation, which describes the unlimited 
(exponential) growth of the scattering amplitude with Y and the symmetric expansion of the support of 
T(Y,r) in r towards both small and large dipole sizes. Note however that the transverse non-locality 
in Eq. (|2.ip is quite weak — the large daughter dipoles with size z ^> r are suppressed by the 'dipole 
kernel' r 2 /[z 2 (r — z) 2 }, whereas the very small ones, with z <C r, are disfavoured by the T-dependcnt 
terms in the r.h.s., which exactly cancel each other as z — > — and can be described as diffusion in the 
logarithmic variable p = ln(rp/r 2 ). Here, rg is an arbitrary scale of reference (say, the unitarization scale 
in the target at Y — 0). 

However, the fully non-linear equation (|2.1|) preserves (and actually saturates) the unitarity bound 
T < 1, as it manifestly has T — 1 as a fixed point. Because of that, the respective evolution is asymmetric 
in r (or p), and the solution T(Y,r) = T{Y : p) looks like a front, which interpolates between T = 1 at 
relatively small p and T = at p — > oo, and which with increasing Y propagates towards larger values of 
p. Behind the front, the scattering amplitude has reached the 'black disk' limit T — 1 and thus cannot 
grow anymore. Ahead of the front, the amplitude is still weak, T <C 1, so the non-linear term in Eq. (|2.1|) 
is unimportant and the amplitude can grow according to the linear, BFKL, evolution. The position 
p s (Y) = ^ D -{ r 'oQ 2 s(X)) °f t ne front at 'time' Y, i.e. the value of p where T becomes of 0(1), represents the 



1 A more detailed discussion of the CCFM equation and its various rewritings will be given in a longer publication [27 



3 



scale where unitarity corrections become important at rapidity Y and defines the saturation momentum 
Ql(Y). 

The previous discussion already suggests that the progression of the saturation front towards larger 
values of p is driven by the BFKL evolution of the dilute tail at p > PsiX) — the front is 'pulled' by 
its tail. This is an important property as it allows us to determine the position p s (Y) of the front and 
its shape around p s {Y) from studies of the linear, BFKL, evolution alone. This property is highly non- 
trivial, in view of the non-locality of the non-linear equation (|2.1[) . Other non-linear equations which 
are non-local and exhibit saturation are known to develop a pushed front, i.e., a front whose progression 
is driven by the growth and accumulation of 'matter' behind the front [281 ]. 

For the case of the BK equation with fixed coupling and for sufficiently high energy, the pulled- 



front property follows from the identification, made in Ref. 15| , between the asymptotic form of the BK 
equation at high energjd and the FKPP equation (from Fisher Kolmogorov, Petrovsky, and Piscounov) 
of statistical physics, for which this property has been established with mathematical rigoJl. However, 
this identification does not extend to a running coupling, and hence it fails to apply for the real QCD 
problem. 

The running of the coupling is known to have dramatic consequences for the high-energy evolution 
SO, QUI, and in particular for the approach towards saturation: the growth of the coupling with 



decreasing momenta, or increasing dipole sizes, amplifies the contribution of the latter to the evolution, 
which then becomes asymmetric even in the absence of saturation. In fact, the BFKL evolution with 
running coupling is infrared-unstable, in the sense that it requires an infrared cutoff to avoid the blow- 
up of the QCD coupling at k± ~ Aqcd, and then the results of the evolution are strongly sensitive 
to the value of this cutoff — so that the whole procedure has no predictive power (see Fig. [3] below). 
In that scenario (BFKL with running coupling), the growth of the gluon distribution at high momenta 
k± 3> Aqcd is mostly feeded by radiation from the bulk of the distribution at infrared (k± ~ Aqcd) 
momenta. Hence, the whole perturbative framework becomes questionable and, besides, one may expect 
the associated saturation front — as generated after enforcing unitarity — to be of the 'pushed' type. 

Yet, as our explicit numerical solutions will demonstrate, this is actually not the case: the saturation 
front remains of the 'pulled' type even for a running coupling. This is so because the gluon modes 
with < Q S (Y) become inert due to saturation, so the evolution is again driven by the dilute tail at 
high momenta, so like for fixed coupling. In particular, the infrared problem is automatically avoided: 
the saturation scale effectively acts as an infrared cutoff, which becomes 'hard' (Q 2 S (Y) ^> Aq CD ) for 
sufficiently high energy. This opens the way towards realistic studies of the front dynamics within the 
context of the linear evolution, as originally suggested in Refs. [I2II13I] . To that aim, the linear evolution 
equations must be supplemented with an appropriate saturation boundary condition, that we now describe. 

Such a boundary condition must ensure that the amplitude never becomes bigger than one. By itself, 
the position of the front is not a priori known, but must be determined when solving the equation. To 
that aim, let us first introduce a line of constant amplitude p = p c (Y) via the condition 

T(Y,p = Pc (Y)) = c, (2.2) 

where the number c is strictly smaller than one, but not much smaller. (The saturation line p s (Y) would 
correspond to c~ 1.) For p < p c (Y) and sufficiently high energy, the solution Tbfkl(Y, p) to the BFKL 
equation would become larger than one — in fact, arbitrarily large. If this equation is to be solved 
numerically, one may think about identifying the point p c {Y) numerically at each step in Y, and then 
enforcing the unitarity limit T = 1 for any p which is smaller than p c (Y) and sufficiently far away from 
it — say, for p < p c (Y) — A with A ~ ln(l/c) a number of 0(1). However, this would not be a very 
good strategy in practice, since T = 1 is not a fixed point for the BFKL equation, so an amplitude of 
O(l) would be exponentially amplified by the subsequent evolution. Even if, at small p, one cuts off this 



2 This is obtained via the gradient expansion of the non-locality in Eq. 112.11 1 to second order in d/dp ('diffusion approx- 
imation'). 

3 The FKPP equation describes a 'reaction-diffusion process' in the mean field approximation corresponding to very large 
occupation numbers at saturation; see e.g. the review paper [28l . 
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evolution by hand step-by-step, it is not clear (especially for running coupling) whether the spurious 
radiation from small p will not affect the tail of the front at large p. It is therefore preferable, as originally 
suggested in Ref. [13] , to enforce the amplitude to vanish for p < p c (Y) — A : 

T(Y,p)=Q for p< Pc (Y)-A. (2.3) 

T = is a fixed point for the BFKL equation, so no further evolution is possible in the 'saturated 
domain' on the left of p c (Y) — A, as it should. When decreasing p below p c (Y), the solution T(Y,p) will 
typically start by rising, then reach a maximum T max ~ 0(1), and eventually decrease to zero. We shall 
conventionally identify the saturation scale p s {Y) with the position of this maximum. In this procedure, 
the numbers c and A are to be viewed as free parameters, which are correlated with each other, since 
A~ln(l/ C ). 

Since it yields, by construction, T — beyond the saturation front, this procedure cannot be used 
for any physical problem which is sensitive to the black disk limit, like deep inelastic scattering at low 
Q 2 < Ql{Y), or particle production at low transverse momenta. On the other hand, as we shall see, this 
procedure accurately describes the dynamics of the front, meaning its position and shape for p > p c {Y), 
and hence it correctly provides the tail of the gluon distribution at transverse momenta k± > Q C (Y), 
with Pc(Y) = \n(r 2 Q 2 (Y)). Since Q C (Y) ~ Q S (Y), we see that our results cover the phenomenologically 
interesting region where geometric scaling is expected at high Q 2 (cf. Introduction). 

To describe our numerical results, let us first change from the coordinate to the momentum repre- 
sentation, i.e. from the scattering amplitude T(Y, r) to the 'unintegrated gluon distribution' A(Y, k) — 
the quantity which enters the calculation of cross-sections within the fey-factorization. For the present 
purposes, A(Y, k) can be defined as the following Fourier transform of the dipole amplitude [ij} 

A(Y,k) = J f^ e - ik - r T(Y,r). (2.4) 

With this definition, the more standard, 'integrated', gluon distribution is obtained as 

AN 2 f® 2 c\ 2 k f 

xg(x, Q 2 ) = j d 2 b A(Y, k, b) . (2.5) 

TT 2 a s J (2tt) 2 J 

For our homogeneous target, A(Y,k,b) = A(Y,k), hence the above integral over b simply yields the 
hadron transverse area irR 2 . With these conventions, the gluon occupation number — i.e., the number 
of gluons of a given color per unit rapidity per unit volume in transverse phase-space — is not exactly 
A(Y, k), but rather A(Y, k)/a s (up to a numerical factor). 

Via Eq. (|2.4p . the saturation front for T(Y,r) translates into a corresponding front for A(Y,k). For 
very large momenta, k ^> Q S (Y), one can use T ~ r 2 ln(l/rf) ('color transparency') and then Eq. (|2.4p 
reproduces the bremsstrahlung spectrum: A(Y, k) ~ 1/fc 2 , ojj A(Y, p) ~ e - ''. This behaviour is modified 
when k gets closer to (but still larger than) Q s , due to the BFKL evolution in the presence of saturation. 
For instance, in the fixed coupling case and for sufficiently high energy (a s Y ^> 1), one finds 12, 13|, 15| 



A(Y,p) ~ (p-pJe-^P-P.) e*p{- (P 2/3a y )2 } fOT Kp-p s <2p s . (2.6) 

The exponential ~ e~ lsP describes the modification of the bremsstrahlung spectrum due to the BFKL 
'anomalous dimension' 1 — j s , while the last, Gaussian, factor describes BFKL diffusion. The numbers 
7s » 0.63 and f3 as 48.52 are specific for the problem at hand: they characterize the BFKL evolution along 
a line of constant amplitude, cf. Eq. (|2.2[) (see 13, 13 1 for details). Also, the presence of the saturation 



scale p s and the overall factor p — p s in the r.h.s. of Eq. (|2.6[) arc the hallmarks of saturation, and can 



4 Prom now on, we shall use the notation p for either ln(r^/r 2 ), or ln(fe 2 /kg) (with ko = 1/ro), the difference being clear 
from the context. 
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Figure 1: The solid lines are the solutions of the BK equation \2.8\ i while the dashed lines have been obtained by applying 
the absorptive boundary J2.3D on the BFKL equation. The leftmost curve has Y = 20 and Y is increased by 10 units for 
each new curve. For the figure on the left we have chosen c = 0.1 and A = 5.0, while for the figure on the right we have 
c = 0.3 and A = 3.0. 



be generated within the framework of the linear evolution only after imposing the saturation boundary 
condition (|2.3p [1 31 ] . Under the same conditions, the saturation scale is obtained as 

pJY) = Xa s Y — hxY + const. (2.7) 

where A sw 4.88 and the last, constant, term is not under control (since Eq. (12. 7|) is merely an asymptotic 
expansion at high-energy). 

In the opposite limit of very low momenta k <C Q S (Y), the integral in Eq. (|2.4p is dominated by large 
dipole sizes r > l/Q s for which T = 1; one then finds A(Y,k) ~ \n[Q s (Y)/k] = [p s (Y) - p]/2. Thus, 
behind the front, A(Y, k) is not exactly constant (unlike the dipole amplitude), but it is slowly growing — 
logarithmically in both 1/k and Y. What saturates at high density is not the gluon occupation number 
A/a s , but rather the rate for gluon emission Q. 

The analytic results in Eqs. (|2.6[) and (|2.7p . together with the corresponding ones at running coupling 
0EI1 j have been already tested in the literature against numerical solutions to the BK equation. Here, 
we are rather interested to compare the solutions to the latter against numerical solutions for the BFKL 
equation supplemented with the boundary condition (12. 3j) — with both fixed and running coupling. The 
momentum-space version of the BK equation is particularly simple in that the non-linear term is local: 

^A(Y,k)=a s J ^ g2(fc 1 _ g)2 (q 2 A(Y,q)- jA(Y,k)j - a s (A(Y,k)) 2 . (2.8) 

The linear version of this equation, i.e., the BFKL equation for A(Y, k), will be solved with an absorptive 
boundary condition similar to that for T(Y, r) in Eq. (|2.3p . This is appropriate since, with the normaliza- 
tion in Eq. (|2.4j) . the saturation effects in the gluon distribution become important when A(Y, k) ~ 0(1). 

For the fixed coupling calculations we shall use a s = 0.2. To include running coupling effects, we 
shall pull the a s factor inside the g-integral in Eq. (|2.8|) and use the one-loop expression for the running 
coupling with scale Q 2 — max(fc 2 ,g 2 ) and Aqcd = 200 MeV. This simple prescription is in agreement 
with the recently constructed running-coupling version of the BK equation 0, E3] ■ To avoid the infrared 
divergence of the coupling at Q 2 = Aq CD , we shall replace a s (Q 2 ) — > a s (Q 2 + p 2 ) for some parameter p. 
Our default choice will be p 2 = 0.5 GeV 2 , but we shall study the sensitivity of our results to variations 
in p. Our initial condition A(Y = 0, k) is given by the bremsstrahlung spectrum for k > 1 GeV (with 
maximal height A — 0.5) and it vanishes for k < 1 GeV. 

For the fixed coupling case, our results are displayed in Fig. [T] for five values of the rapidity within 
the range 20 < Y < 70 and for two different choices for the parameters c and A. The most important 




k (GeV) k (GeV) 

Figure 2: The running coupling results for: BK (solid curves), BFKL with absorptive boundary (long dashed curves) and 
pure BFKL (short dashed curves) for (left) Y = 10,20,30 and 40, and (right) Y = 6,8, 10, 12 and 14. For the absorptive 
boundary we used c = 0.1 and A = 5.0. 



observation about these results is that the saturation fronts generated by the two types of evolution do 
precisely coincide with each other for all momenta p > p c (Y) — A, and for all the considered rapidities. 
This property is not altered by changing the values for c and (correlated to it) for A. We have also 
checked that these numerical curves are consistent with the analytic estimates in Eqs. (|2.6p and (|2.7[) . 
including the expected values for A and j s . 

We now turn to the more realistic case of a running coupling. Then, as alluded to before, the pure 
BFKL evolution is infrared unstable, since the rapid growth of the gluon distribution at small values of k 
(where is coupling is larger) is also feeding the growth at higher k. Therefore the linear evolution behaves 
quite differently compared to the non-linear one, even at high k. This is clearly visible in Fig. [5] where 
we compare the strict BFKL evolution to the BK one, and to the BFKL evolution with the absorptive 
boundary. Once again, there is a perfect matching between the saturation fronts provided by BK and, 
respectively, BFKL with saturation boundary. (For the latter, we used the same values for c and A as 
at fixed coupling.) On the other hand, we see a dramatic difference with respect to the linear evolution, 
which progresses much more rapidly towards the right. 

From these curves, it is also possible to extract the K-dependence of the saturation momentum p ? 
for running coupling. We find that the squared-root law p s ~ X r \/Y predicted by the theory [IT 
for asymptotically high energies provides a good fit to our numerical results for Y > 10, with a fitted 
value A r ~ 2.8 which agrees reasonably well with the (asymptotic) theoretical expectations A r ~ 3.2. 

Now, from the phenomenological point of view, we are more interested in values of Y which are not 
that large, say Y < 14 (corresponding to x > 10~ 6 ), as relevant for forward jet production at LHC. With 
that in mind, we also show in Fig. [2] (right) the results for lower values of Y, between 6 and 14 units; one 
can thus see that the absorptive boundary method works equally well also for such lower rapidities. 

Finally, to illustrate the infrared stability introduced by saturation, we exhibit in Fig. [3] results 
obtained for different values of the IR cutoff p 2 inserted in the running coupling. Unlike the pure BFKL 
results (left figure), which are extremely sensitive to a change in p, the results corresponding to the 
saturation boundary condition (right figure) show no sensitivity whatsoever. 



° Specifically, for asymptotically large Y , the running -coupling BFKL evolution yields [H,[l]| : ps(Y) ~ y/2\b Y where 
A ~ 4.88 is the same number as in Eq. <2~7t and b = 127V C /(117V C - 2N f ) is the coefficient in the one-loop running 
coupling: ct(Q 2 ) = bo /ln(Q 2 / Aq CD ) . In our simulations, we use Nf = 0, hence we expect A,- = \/2A&o — 3.26, which is 
indeed consistent with the fit to the curves in Fig. [2] 
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k (GeV) k (GeV) 

Figure 3: Testing the sensitivity of the BFKL evolution to the IR cutoff fi. Left: The pure BFKL evolution with u? = 0.5 
GeV 2 (solid lines), u 2 = 1.0 GeV 2 (long dashed lines) and p? = 2.0 GeV 2 (short dashed lines), for Y = 10, 20, 30 and 40. 
Right: The BFKL evolution with absorptive boundary for u 2 = 0.5 GeV 2 (solid lines), u 2 = 2.0 GeV 2 (dashed lines) and 
for the same values for Y as before. 

3. CCFM evolution with absorptive boundary 

In this section we shall present a compact version of the CCFM equation Q to which we shall apply 
the boundary condition described in the previous section. A more comprehensive discussion of the CCFM 
formalism and its relation to BFKL will be given elsewhere [531, together with more detailed numerical 
studies, of which the present Letter is only giving a glimpse. 

As mentioned in the introduction, the CCFM evolution takes into account the quantum coherence 
between successive emissions via angular ordering in the parton cascades. Accordingly, the respective 
gluon distribution now depends on three variables, A — A(x, k, q), where the third variable q is a trans- 
verse momentum related to the maximum angle which determines the phase space where emissions are 
allowed. This angle is set by the hard scattering of the space-like photon against a quark inside the 
proton. It is customary to define the variable £ which is the squared angle, £ = q 2 / (y 2 E 2 ), where q is the 
transverse momentum of a gluon emitted in the s-channel, y is its longitudinal momentum fraction, and 
E is the energy of the proton; then, all emissions must satisfy £,<£,= q 2 /(x 2 E 2 ). 

The CCFM equation for A can be written in different versions, depending on how 'exclusive' we 
choose the gluon distribution to bqj. That is, so long as one is not interested in the structure of the final 
state, one can 'integrate out' some of the emissions in the s-channel, as they do not change the overall 
(unintegrated) gluon distribution (and hence neither the probability for the interaction with a projectile). 
In practice, this amounts to suitable cancelations between 'real' gluon emissions and 'virtual' terms (or 
'Sudakov factors'). Of course, if one is interested in studying the exclusive final states, all emissions which 
were removed from the initial state must later be included as final state radiation. 

If one keeps within A only the emissions associated with the \ j z pole in the splitting function, then 
the CCFM equation can be written as the following integral equation 

r dz f d 2 n fx \ 

A(x,k,q) = a s — —§0(q-zq)A ns (k,z,q)A{~,\k + q\,p), (3.9) 

where q and 1 — z are, respectively, the transverse momentum and the energy fraction of the 'real' gluon 
emitted (in one step of the evolution) in the s-channeQ- The theta function comes from the angular 
ordering constraint £ < £. There is also an energy ordering implicit in (|3.9|) : emissions are ordered in 



6 Recall that A(x, k) represents the distribution of gluons produced by the evolution in the t-channel, as available for 
the interaction with an external projectile. 

7 Strictly speaking, Eq. 1 13. 91 1 should involve the rescaled variables p = q/(l — x) and p = q/(l — z), but here we are only 
interested in the small-x behaviour, so we make no distinction between e.g. p and q [27H . 
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energy as well as in angle. Finally, A ns is the so-called 'non-Sudakov form factor', which accounts for 
virtual corrections and is necessary to ensure probability conservation 

A ns (z, k,q) = exp (-®s J J ^ ^ = Cxp ( _ " s log Q) log (jf )) - ( 3 ' 10 ) 

Since A now depends upon three variables, Eq. (|3.9p is much more difficult to solve than the BFKL 
equation. To simplify the numerics and thus be able to explore a relatively wide range in k and Y, it 
is convenient to use a simpler version of the CCFM equation, that we shall now derive. This is possible 
since, as first noticed in [25| (see also Ref. [Ill), one can further integrate out some of the real emissions 
explicit in Eq. (|3.9p . and thus get a more inclusive equation. The main point, to be further detailed in 



27] . is that the non-Sudakov factor in (|3.10p , despite its name, can be used as a kind of a Sudakov factor, 



to cancel a certain class of real emissions. 

Since in CCFM (like in BFKL) we integrate over the transverse momentum q of the emitted gluon, 
one can roughly have three possibilities: either k' w k ^> q, or k' w q ^S> k, or, finally, k sw q 3> k 1 . (Here 
k' = \k + q\, so we consider the elementary splitting k' — > k + q.) The emissions satisfying the first 
condition are the ones which can be canceled against A ns . To that aim, we must also include, within 
the integrand of Eq. (|3.9|) , the so-called kinematical constraint k 2 > zq 2 which ensures that the squared 
four-momenta of the virtual propagators are dominated by their transverse part as required by the multi 
Regge kinematics. This constraint is generally kept implicit in the CCFM (or BFKL) literature, since 
the Regge kinematics is guaranteed to the order of interest; yet, its explicit inclusion in the equations 
introduces corrections which are formally of higher order in a s , but which can be numerically important. 
After including this constraint, we can remove the factor A„ s from Eq. (|3.9[) and simultaneously limit 
ourselves to emissions satisfying the last two constraints written above, which can be summarized as 
9(q 2 — min(/c 2 , k 12 )). We thus obtain 

A(x, k,q)=a s I — f^9(q- zq)9(k 2 - zq 2 )9{q 2 - mm(k 2 ,k /2 ))A (—, k', q) . 

(3.11) 

Since q > k for all cases of physical interest, we further have q 2 > k 2 > zq 2 > z 2 q 2 . Therefore the 
angular ordering is automatic and 9{q — zq) can be neglected. This means that the dependence on the 
third variable q drops out, and we can write (note also that q > k') 

A(x, k) = a s f 1 — 9(k 2 - zq 2 )9(q 2 - mm(k 2 ,k' 2 ))A (-, k') , (3.12) 

Jx z J \z J 

which as compared with the original Eq. (|3.9p represents a considerable simplification. 

We shall now perform the integration over the azimuthal angle <j> between q and k. To that aim, it 
is convenient to replace 9(k 2 — zq 2 ) by 9(k 2 — zk' 2 ), which is allowed within the current approximation^, 
and then switch the integration variables from q to k' and, respectively, from z to x/ z (which we rename 
as z). This yields 

A(x,k)=a s f 1 — P ^ f ,J k ' 2 ., 9 9(z - xk' 2 /k 2 ) 0(\k' - k\ 2 - min(fc 2 , k' 2 ))A(z, k'). 
Jx z Jo 77 J \k — k\ 



(3.13) 



Notice that this equation is infrared finite. After also performing the 4> integral, one finds 



r 1 c\? r t\k' 2 

A(x, k)=a s I — 9{z - xk' 2 /k 2 ) h( K ) A(z, k') , (3.14) 



8 With this we mean that the total probability to have any number of real emissions of the given class times the form 
factor is one (the emissions are probability conserving). 

^Indeed, we have either k' q S> k, in which case the replacement is obviously correct, or k m q 2> k' , in which case 
both the first and the second theta function can be replaced by 1. 
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k (GeV) k (GeV) 

Figure 4: Comparison between the solutions to the CCFM equation l|3.16ll (solid lines) and the BFKL equation (dashed 
lines) for Y = 10, 20, and 30. Left: fixed coupling a s = 0.2. Right: running coupling. 




where k = min(fe , k )/max(fc , k ) and 



W . 1 - i«rct,„ 9( « - 1/4 ). (3.15) 

Notice that h(n) — > as k — > 1, so Eq. (|3.14[) has indeed no singularity at fc' = fc. To obtain an integro- 
difFerential equation, we define z = e~ v and x — e~ Y and differentiate the l.h.s. with respect to Y. We 
thus get 

r c\k' 2 / 
d Y A(Y,k) = a s J |fc2 _ h(K) \9(k 2 - k' 2 )A(Y, k') 

+6(k' 2 - k 2 )9{Y - \og(k ,2 /k 2 ))A(Y - log{k' 2 /k 2 ), k')) , (3.16) 



which is our final version for the CCFM equation and is equivalent to an equation originally proposed in 
Ref. [5^1 (although there this was not derived from the CCFM equation (|3.9[) in the way we did). The 
fact that this equation is nonlocal in Y should not come as a surprise, since the CCFM evolution is not 
really an evolution in Y , but rather in £. 

It is first interesting to compare the predictions of Eq. (|3.16[) to those of the BFKL equation for the 
strictly linear evolution. This is shown in Fig. 3] for both fixed and running coupling, and one clearly sees 
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Figure 6: Left; Solutions to the running-coupling CCFM equation H3.16J I with the saturation boundary (solid lines) and 
without it (dashed lines) for Y = 8, 10, 12, and 14. Right: CCFM (solid lines) vs. BFKL (dashed lines) solutions with 
running coupling and saturation boundary for very high rapidities: Y = 60, 70, 80, 90, 100, and 120. 



that the BFKL evolution is considerably faster. This difference is to be attributed to the non-local term 
in the r.h.s. of Eq. (|3.16|) : the 'retarded' distribution A(Y — log(fc' 2 /k 2 ), k') is generally smaller than the 
'instantaneous' one A(Y,k'). 

But even though the CCFM evolution is somewhat slower, Eq. (|3.16p still shows a pronounced 
growth with Y, which in the absence of any non-linearity would rapidly lead to unitarity violation. To 
cure for that, we now enforce the absorptive boundary condition (|2.3p on Eq. (|3.16|) . The corresponding 
results are compared to those of the purely (CCFM) linear evolution in Fig. [5j for both fixed and running 
coupling. As in the BFKL case (compare to Fig. [2|), the difference is more pronounced for a running 
coupling, since then the linear evolution is again infrared-unstable, and this instability is removed by the 
inclusion of saturation. 

For the more realistic, running-coupling, case it is furthermore interesting to show the results at 
lower rapidities, as relevant for LHC. This is exhibited in the leftmost figure in Fig. [51 together with the 
corresponding results of the strictly linear evolution. As one can see there, for Y = 14 and k as high as 
10 GeV (which is well above the respective saturation momentum Q s ~ 2.5 GeV), saturation reduces the 
predicted gluon distribution by about a factor of 2. 

Finally, in the rightmost figure in Fig. [6l we compare the saturation fronts provided by the BFKL 
and CCFM evolutions with running coupling and for relatively high rapidities, up to Y = 120. What is 
remarkable about this last figure is that the BFKL and CCFM evolutions with saturation and running 
coupling appear to very close to each other and for all rapidities — meaning that the corresponding 
fronts propagate at roughly the same speed. This is at variance with the corresponding situation at fixed 
coupling, where the BFKL evolution is still faster. This is confirmed by the estimates of the saturation 
momentum for the CCFM evolution, as extracted from fits to our numerical results: It is numerically clear 
that p s rises linearly with Y for a fixed coupling, and it exhibits a VY behaviour for a running coupling, 
just like for BFKL. Trying simple fits of the form p s = Xfa s Y and, respectively, p s — X r VY, we find the 
value J^l A / ~ 3.5 and, respectively, A r ~ 2.8. For fixed coupling, this value A/ for the saturation exponent 
is indeed smaller (although not much smaller) than the corresponding BFKL estimate A/ ~ 4.9. But for 
running coupling, the CCFM and BFKL estimates for A r are essentially the same within our numerical 
accuracy. Moreover, the manifest similarity between the shapes of the BFKL and CCFM fronts in Fig. [6] 
(right) suggests that the CCFM fronts exhibit the same properties of geometric scaling as the BFKL ones 
(again, for a running coupling). 



10 Note that, at fixed coupling, the CCFM estimate for the parameter Xf defined as above is still a function of a s (unlike 
the respective BFKL estimate); hence this value Aj ss 3.5 must be seen as the value corresponding to a s = 0.2. We shall 
further discuss the o s -dependence of the CCFM parameter Aj in Ref. [27l . 
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It would be of course interesting to understand this similarity between the BFKL and CCFM evolu- 
tions towards saturation in more depth, and also to perform more detailed studies of the CCFM evolution 
with saturation boundary, in order e.g. to explicitly test geometric scaling. We postpone such studies to 



a further work [27 1 
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